/*
 * AqueousPhase.h
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#ifndef AQUEOUSSOLUTION_H_
#define AQUEOUSSOLUTION_H_

// Eigen includes
#include <Eigen/Dense>
using namespace Eigen;

// GeoReact includes
#include "AqueousActivityModel.h"
#include "Utils.h"

class AqueousPhase
{
public:
	enum IonicActivityModel { Davies, ExtendedDebyeHuckel, Hkf, Pitzer };
	
public:
	AqueousPhase();
	
	explicit AqueousPhase(const vector<string>& species);
	
	void
	SetIonicActivityModel(IonicActivityModel model);
	
	void
	SetSetschenowConstant(const string& speciesName, const vector<double>& bi, const vector<double>& Ti);
	
	const vector<string>&
	GetSpecies() const;
	
	const vector<string>
	GetIonicSpecies() const;
	
	const vector<double>
	GetIonicCharges() const;
	
	const Index
	GetWaterIndex() const;

	const double
	EffectiveIonicStrength(const vector<double>& n) const;
	
	const double
	StoichiometricIonicStrength(const vector<double>& n) const;
	
	const double
	TotalMolalitySolutes(const vector<double>& n) const;
	
	const vector<double>
	Molalities(const vector<double>& n) const;
	
	const vector<double>
	IonicMolalities(const vector<double>& n) const;
	
	const vector<double>
	IonicStoichiometricMolalities(const vector<double>& n) const;
	
	const vector<double>
	ActivitiyCoefficients(double T, double P, const vector<double>& n) const;
	
	const vector<double>
	Activities(double T, double P, const vector<double>& n) const;
	
private:
	/// The species in the aqueous phase.
	vector<string> species;
	
	/// The index of the water species.
	Index iwater;
	
	/// The vector of electrical charges of the ionic species.
	vector<double> zi;
	
	/// The indexes of the ions, neutral species and electrolytes.
	vector<Index> ions, neutrals, electrolytes;
	
	/// The stoichiometric matrix between the electrolytes and ions.
	MatrixXd vei;
	
	/// The function that computes the activity coefficients of the ionic species.
	IonicActivityFunction ionicActivityFunc;
	
	/// The function that computes the activity coefficients of the neutral species.
	vector<NeutralActivityFunction> neutralActivityFuncs;
	
private:
	void 
	InitElectrolytes();
	
	void 
	InitDefaultActivityModels();
	
	const vector<double>
	ActivityCoefficientsIonicSpecies(double T, double P, const vector<double>& n) const;
	
	const vector<double>
	ActivityCoefficientsNeutralSpecies(double T, double P, const vector<double>& n) const;
};

#endif /* AQUEOUSSOLUTION_H_ */

